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ABSTRACT 

An artificial dissipation model, including boundary treatment, that is 
employed in many central difference schemes for solving the Euler and Navier- 
Stokes equations is discussed. Modifications of this model such as the eigen- 
value scaling suggested by upwind differencing are examined. Multistage time 
stepping schemes with and without a multigrid method are used to investigate 
the effects of changes In the dissipation model on accuracy and convergence. 
Improved accuracy for Invlscld and viscous airfoil flows is obtained with the 
modified eigenvalue scaling. Slower convergence rates are experienced with 
the multigrid method using such scaling. The rate of convergence is Improved 
by applying a dissipation scaling function that depends on mesh cell aspect 
ratio. 
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!• INTRODUCTION 


In the past few years, substantial progress has been achieved in the 
development of efficient numerical schemes for solving the Euler and Navier- 
Stokes equations^""^* Robustness and accuracy of the schemes has also con- 
tinued to improve. Strong emphasis has been placed on sharp representation of 
shock waves, which is reflected in the Euler solutions obtained^""^®. Now, the 
accuracy of viscous flow calculations . (i.e. , turbulent flows where there are 
strong gradients) requires additional attention. For example, nonphysical 
solutions have been obtained for trailing edge turbulent airfoil flows^^"^^. 
A major factor contributing to inaccuracies is the artificial dissipation 
present in the numerical algorithms. 

The schemes that are used for solving the Euler and Navier-Stokes equa- 
tions are based on either central or upwind differencing. Both central and 
upwind methods include artificial dissipation. A S3rmmetric form^^ for the 
numerical flux function clearly reveals that upwind schemes involve a matrix 
dissipation coefficient. This results in a specific scaling (based on charac- 
teristic values) of the dissipation of each conservation equation. In the 
case of central difference schemes, a scalar coefficient is employed for the 
dissipative flux contribution to the numerical flux. This results in a 
simpler scheme with a smaller operations count. For either type of dif- 
ferencing, the principal requirements in the design of the dissipative terms 
are that they must be large enough for a satisfactory convergence rate and yet 
sufficiently small that accuracy is not compromised. 

In this paper, a central differencing algorithm is used to investigate 
artificial dissipation. There are two fundamental reasons for adding dissipa- 
tion terms to a central difference method. First, they are included to pro- 
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vide high frequency damping. It is well-known that central difference schemes 
experience odd and even point decoupling for both linear and nonlinear prob- 
lems. These high frequency modes must be damped to achieve satisfactory 
convergence. In the case of nonlinear problems, high frequency damping is 
required to remove the energy produced by nonlinear interactions (i.e., con- 
sider a Fourier representation of nonlinear convection terras). Without such 
damping, the unresolvable modes (subgrid frequency components) can appear as 
errors in the resolvable low frequency components of the discrete solution. 
Second, artificial dissipation terms are added to eliminate oscillations in 
the neighborhood of shock waves. Also, from the mathematical theory for 
hyperbolic systems of inviscid conservation laws^^, the introduction of 
artificial dissipation is necessary to guarantee a unique weak solution. 

It is interesting to note that if sufficient resolution were used to 
define a shock structure, the solution of the full Navier-Stokes equations 
would eliminate the need for artificial dissipation at shock waves^^. How- 
ever, this would mean that the mesh spacing in the streamwise direction in the 
vicinity of the shock would have to be orders of magnitude (depending on the 
Reynolds number) smaller than that which is currently used in aerodynamic 
computations. Furthermore, solving the complete Navier-Stokes equations 
rather than a subset such as the thin-layer Navier-Stokes equations (where 
diffusion terras in the streamwise-like direction are neglected) could require 
much greater computer time. 

In the present work, the artificial dissipation model Introduced by 
Jameson, Schmidt, and Turkel^^ is reviewed. Then, some modifications of this 
model and boundary treatment of the dissipative terms are discussed. Numeri- 
cal methods used to solve the Euler and thin-layer Navier-Stokes equations are 
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briefly described* Next, inviscid, laminar, and turbulent airfoil flows are 
considered to investigate the effects of certain modifications of the basic 
dissipation model on efficiency and accuracy. Special emphasis is given to 
the calculation of accurate viscous flow solutions. 


II. BASIC DISSIPATION MODEL 

The basic dissipation model considered in this paper was first introduced 
by Jameson, Schmidt and Turkel^® in conjunction with Runge-Kutta explicit 
schemes. It has subsequently been used by many lnvestlgators^^“^^ in a wide 
range of applications. Also, it has been applied to ADI implicit schemes^^. 
In this section, this model will be briefly reviewed. 

Consider the Euler equations in the form 

Wt + fx + gy = 0 (1) 

where W is the solution vector of conserved variables, and f, g are the 
inviscid flux vectors. The Independent variables are time t and Cartesian 
coordinates (x,y). Transforming Eq. (1) to arbitrary curvilinear coordi- 
nates C = 5(x,y), n = n(x,y) 

(J • W)^ + = 0 (2) 

where J is the transformation Jacobian and F = fy^ - ^ * 

In a cell centered finite-volume method, Eq. (2) is simply integrated over an 
elemental volume in the discretized computational domain, and J is then the 
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volume of the cell. Eq. (2) can also be written as 


+ AW_ + BW =0 
t 5 n 


where A and B are the flux Jacobian matrices. 

A typical step of a Runge-Kutta approximation to Eq. (2) is 


- a, j”^At [D + D - D] 

K S T1 


where D^, are approximations to the spatial derivatives, and D are 
artificial dissipation terms, which are usually frozen at the first or second 
stage. The artificial dissipation employed in Ref. 18 is a blending of second 
and fourth differences. That is. 


D = (Dp + - D^)W 

? n 5 n 


)^W • ^i+l/2,j^^C^i,j 




and A^ , 7^ are forward and backward difference operators associated with 
the 5 direction. The variable scaling factor 


^1+1/2, .i ^^C^i+l,j 


where 


is the largest eigenvalue of the matrix A, and 


is the 
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largesC eigenvalue of the matrix B. The coefficients and 

are adapted to the flow and are defined as follows: 



where P is the pressure, and typical values of the constants and 

are 1/4 and 1/256, respectively. The operators in Eq. (5) for the 
n direction are defined in a similar manner. 

Before proceeding, some general comments on the form of these terms are 
appropriate. First, the original use of this artificial dissipation was for 
the solution of the Euler equations on a grid with an aspect ratio close to 
one. Second, the scaling factor X, which is given in Eq. (8), has an iso- 
tropic behavior. Such a behavior is generally not satisfactory in viscous 
flow calculations. Also, the eigenvalues 



where u, v are Cartesian velocity components and c is the speed of sound, 
represent approximations to the flux Jacobian matrices A and B. (See Refs. 
15 and 24 for the relationship between central differencing plus artificial 


i 
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dissipation and upwind differencing.) Finally, in more recent versions of the 
dissipation model, the maximum in Eq. (9) is taken over more mesh cells than 
the immediate neighbors. This is beneficial for shock capturing capability. 


III. MODIFICATIONS OF BASIC DISSIPATION MODEL 

The second difference dissipation term given in Eq. (6) is an approxima- 
tion to 

k lf>- 


where = X • Adding this expression to the right-hand side (RHS) 
of Eq. (2), multiplying the resulting equation by W, and integrating over the 
domain (n) gives 


1/2 • Jdgdn = flux terms 

(13) 

r ( Al) ,3W,2 


if boundary terms are neglected or if boundary derivatives vanish. For linear 
problems, the square of the norm j j • Jd^dn (which in this case 

a 

is an energy estimate in the mathematical sense) is a good measure for the 
stability of the numerical scheme. Equation (13) shows that the second dif— 
ference dissipation term decreases this norm and, thus, is strictly dis- 

sipative. If the same type of analysis is done for the fourth difference 
dissipation term of Eq. (7), then 
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1/2 J / • Jd^dn = flux terms 

-If 04 ) 

ft ac 

f f ag^^^aw. a^w ... 

ft as ^ as 

Both a dissipative term and a dispersive term appear on the RHS of Eq. (14). 
The following terra 


■>«" - (h,i 


is considered as a replacement for the one in Eq. (7). This modified term 

(4) 

produces only dissipative terms. Note that X and e are evaluated 

at nodes rather than at mesh cell boundaries as in Eq. (7). 

For Navier-Stokes problems, a fine mesh is required in the direction 

normal to the body in order to resolve the boundary layer. In the interest of 

computational efficiency, the mesh spacing in the streamwise direction for 

high Reynolds number calculations is generally chosen so as to resolve the 

streamwise inviscid terms only (i.e., thin-layer Navier-Stokes assumption). 

Then the mesh in the viscous region has a high aspect ratio (with S as arc 

length AS /^S^ « 1). To make matters more difficult, the situation can be 

reversed in the far field of an external flow problem. Thus, depending on the 

grid generation technique AS /AS^ =0(1) or even AS /AS^ » 1 in the 

n 5 n 5 

far-field region. These large distortions create difficulties both for the 
convergence and for the accuracy of steady-state computations. These 
difficulties are compounded for multigrid schemes since high frequency modes 
are very different in the two coordinate directions. 


J 
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A number of investigators have suggested that an anisotropic dissipation 
model is needed for such problems. Therefore, Eq. (8) in the basic dissipa- 
tion model is replaced by 


= 1/2[(X_). . + 
^ ^ » J 


and a similar equation is used in the n direction. For a multigrid 
algorithm, this scaling in the streamwise direction can be too severe. More- 
over, the effectiveness of the driving scheme in damping high frequencies in 
the C direction can be significantly diminished, resulting in a much 
slower convergence rate. In Ref. 25, Martinelli introduces functions of mesh 
cell aspect ratio and obtains accurate solutions and good convergence rates. 
For example, one can replace Eq. (16) by 



optimum 
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Due to large velocity gradients in turbulent boundary layers, additional 
scaling of the artificial dissipation is required in the direction normal to a 
surface boundary. The presence of the physical viscous terms can be exploited 
to allow the additional dissipation terms in the normal direction to be 
reduced. In the present work, this is accomplished by multiplying the second 
and fourth difference dissipation terms by a simple linear function of the 
local Mach number. That is, the normal scaling factor becomes 


^l,j+l/2 



n i,j 


^\^l,j+l^ 


( 18 ) 


where 


M = and 


is the local Mach number. 



IV. BOUNDART TREATMENT OF DISSIPATION TERMS 

In this section, the boundary dissipation operators that are applied in 
many flow prediction codes based on finite-volume discretization are pre- 
sented. Then, a local mode analysis is used to examine the relative damping 
characteristics of some of the difference stencils. The Influence of the 
boundary cell operators on the character of the dissipation matrix for the 
system of flow difference equations is also discussed. 

In a finite-volume method, the first and last cells in each coordinate 
direction are auxiliary cells where the flow equations are not solved. The 
solution in these cells is found by a combination of the given physical 
boundary conditions and numerical boundary conditions (l.e., extrapolation). 
Hence, there is generally no difficulty in evaluating the second difference 
dissipation term at the first or last interior cell in a given coordinate 
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direction# Note that at a solid surface boundary either the surface or entire 
contribution to the normal-like dissipation of the first interior cell is 
usually set to zero# In the case of the fourth difference dissipation term, 
information is required at two neighboring cells on each side of the cell 
being considered. Therefore, special treatment of this term is needed for the 
first interior cell at the boundaries of the physical domain. Eriksson and 
Rizzi^^ and Pulliam^^ suggest choosing a boundary cell difference stencil that 
results in a nonpositive definite dissipation matrix for the system of dif- 
ference equations. As will be shown, such a choice results in a numerical 
scheme that is more dissipative for the long wavelength components of the 
solution at a boundary than in the interior of the domain. Although this may 
be acceptable at a far-field boundary of an external flow problem, caution 
should be exercised in selecting the difference formula at a solid boundary. 
For example, in Euler calculations a large dissipation in the direction normal 
to the boundary can generate a thick false entropy layer. Also, as indicated 
previously, it can alter a viscous flow solution significantly. 

At this point some simplifying notation is introduced to identify and 
subsequently analyze some of the boundary cell treatments that are commonly 
used for the fourth-difference dissipation. First, let D, F, F, and G 
denote first, second, third, and fourth differences, respectively. Then, for 
interior cells in a given direction 

Vl/2 ■ \ 

h - Vl/2 - Vl/2 - Vl - ^ Vl 
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^+ 1/2 ■ - '=11 ■ V 2 - - “ n -1 


‘=2 - ^+1/2 - ^-1/2 ■ V2 - ^ % - %-l ^ V2- 


where is the discrete solution for the cell. The dissipation 
stencils considered for the first interior cell (designated 1=2) are 
generated by applying the following (see Fig. (1)): 


(Al) ^1/2 ~ ^3/2 (zeroth order extrapolation) or = 0. 

(A2) ^1^2 ~ ^^3/2 ~ ^5/2 (first order extrapolation) or = E 2 > 

(A3) = °3/2 °5/2 " °7/2 " ^3* 

(A4) = ^^^3/2 ~ ^5/2^ ^7/2 (quadratic interpolation) or 

Ef = -F ^^2 + °5/2 “ °3/2* 


In the case of a solid surface boundary (l = 3/2), the normal difference 
operators that are generally used are constructed by setting the surface dis- 
sipative flux ^ 3/2 zero and 


<B 1 ) D 3/2 ■ “ 5/2 ®2 - ”• 

(B2) - 2Djy2 - ^2 ' ®3' 

(B 3 ) ^ 5/2 “ 0 > ^3 “ 0 (numerical dissipation of zero for first two 

interior cells). 


The treatment of (B3) has been applied successfully in both inviscid and 

11 23 

viscous multidimensional flow calculations. ^ 


- 12 - 


A local mode analysis can be beneficial In examining the relative damping 

behavior of boundary cell difference operators. First, for comparison 

purposes, we characterize the interior fourth difference. Taking the Fourier 

transform of G. we obtain 
I 

g^ = 4(cose-l)^ 

where g^ is the Fourier symbol and 0 is the product of a wave number 
and the mesh spacing. Then, 


— 4 

g^ ~ 0 for small 0 

and 

g^Cir) = 16. 

The dissipation of long waves is dictated by the behavior of at small 

0, and the dissipation of short waves is governed by The 

X 

coefficient (see Eq. (11)) is chosen so that the highest frequency is 

highly damped. This is important for multigrid calculations. Near a 
boundary, the dissipation should behave in a similar manner. 

A general form of the difference stencil at i - 2 and the associated 
Fourier transform symbol can be written as follows: 

\ - * (6 ♦ If - a)W^ - ifW,., t . 2 

g^ = [6 + Y - 2a(l + cos0)](l - cos0) 


and 
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+ 1(y - 6 + 2a cos9)sln9, 
gjj “ (3+Y-4a)j + o|- 

+ i(2a - 6 + y) 0 ” io9^ for small 9 
g^(Tr) = 2(3 + y). 


In the case of (Al) 


‘ ■ 2 . 


16, 


g. “ 4(1 - cos9) - 2e (1 - cos9). 


— 23 

g^ « 9 - 19 for small 9 , 


gj^(lT) = 12. 


Note that is not real; and thus, there is both dissipation and 

dispersion near the boundary. This is the dissipation recommended by 
Pulliam^^. It is second order on long waves and fairly dissipative on short 
waves. The treatment of (A2) gives 


=1 ■ - 3Vi ^ - “i-i * * 


=2(1 - COS0) - 2isln6(l - cos0), 
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. 0 2 

® 2 small 9 , 

= 8 . 

For all waves, the Real (g. ) is half of the interior value. The dtssi- 

X 

pation formula for (A3) is simply twice that of (A2), Then, the real part 
of for all waves becomes the same as it is in the interior of the do- 

main. In numerical experiments, the boundary cell treatments of (Al) - (A4) 
resulted in similar solutions and convergence rates. 

In Ref. 24, boundary difference stencils are evaluated by computing the 
eigenvalues of the dissipation matrix for a one-dimensional discrete system 
that includes a fourth difference dissipation term. Such an evaluation shows 
that the damping of the highest frequency (as determined by the largest eigen- 
value, say A ) by boundary treatments (Al) - (A3) is nearly the same, 
max 

Note that as the number of mesh points increases , the eigenvalue X is 

max 

dictated by the interior point stencil (for the interior “ 16). The 

principal difference between using (Al) or (A2) - (A3) is indicated by the low 
frequency behavior of the dissipation matrix (see also small 

0). The matrix associated with (A2) or (A3) has a zero eigenvalue. There- 
fore, (A2) or (A3) are not recommended since they could lead to undamped 
modes. According to Refs. 26 and 24 a boundary dissipation formula is chosen 
so that the dissipation matrix is nonpositive definite (i.e., strictly dissi- 
pative). The dissipation treatments resulting from (Al) and (Bl) satisfy this 
requirement. The results of this paper were obtained with (Al) and (Bl). 
Moreover, the boundary treatment and interior representation (Eq. (15)) of the 
fourth difference dissipation are consistent (i.e., strictly dissipative). 
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V. NUMERICAL METHODS 

The numerical results presented in this paper were computed with multi- 
stage time stepping schemes. Details and properties of these schemes have 
been described previously^ ^ , In some of the calculations, both four and five 
stage Runge-Kutta algorithms were used as drivers for a multigrid process. 
The multigrid technique is based on the work of Jameson. In particular, a 
Full Approximation Storage (FAS) method^® and V-type cycle are employed. The 
grid transfer operators (i.e., restriction and prolongation operators) are the 
same ones used in the Jameson procedure. Several modifications of the 
original method have resulted in improved multigrid performance. First, the 
fourth difference dissipation terra is coraputed with Eq. (15). The norraal 
artificial dissipation near the wake line is treated by continuation rather 
than applying the sarae procedure used on the airfoil surface. All boundary 
inforraatlon is updated after each stage and on all meshes in the multigrid 
process. Finally, on each level of refinement of a Full Multigrid (FMG) 
method, multiple iterations are performed on coarse grids. One iteration is 
done on the finest mesh, two Runge-Kutta cycles on the next mesh, and three 
Runge-Kutta cycles on all coarser meshes. In the viscous flow calculations, a 

OQ 

convective coarse grid correction scheme is used. Moreover, the viscous 
terms are evaluated only on the finest grid for a given level of refinement. 
For further discussion of the multigrid algorithm, see Ref. 30. 


VI. RESULTS AND DISCUSSION 

Adequate consideration must be given to convergence as well as accuracy 
in designing an artificial dissipation model. This is especially true for a 
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tnultigrid technique. For example, good high frequency damping of the basic 
solver is the crucial requirement for constructing an efficient multigrid pro- 
cess. In the first part of this section, the effects of scaling of the 
numerical dissipation are investigated by considering multigrid calculations 
for inviscid and laminar flow over an airfoil. The last part deals with tran- 
sonic turbulent airfoil flow, and in particular, the trailing edge flow. To 
facilitate the discussion of the numerical results, the following designations 
are made to indicate the form of the artificial dissipation model used: 

1) Basic or original (see Section II) 

2) Modified (Eq. (16)) - This refers to scaling with individual 

eigenvalues. 

3) Modified (Eq. (17)) - Individual eigenvalues are multiplied by a 

function of cell aspect ratio. 

For each model, Eq. (15) is used for the fourth difference dissipation term. 
Finally, wherever a convergence history is presented, it shows the variation 
of the logarithm of the root mean square of the residual of the continuity 
equation with Iteration. For the multigrid computations, an Iteration 
corresponds to a multigrid cycle. 

Transonic Inviscid Flow 

Several calculations were performed for an NACA 0012 airfoil at Mach 0.8 
and an angle of attack of 1.25°. A C-type mesh with 256 cells around the air- 
foil (193 points on the airfoil) and 32 cells normal to the airfoil was 
used. The outer boundary was placed 12 chords away from the airfoil; a far 


-17- 


O t 

field vortex boundary condition'^^ was applied. The computed surface pressure 
distributions and convergence histories using the basic and modified (Eq. 
(16)) artificial dissipation models are displayed in Figs. 2 and 3, respec- 
tively. The predicted shocks on the upper and lower surfaces of the airfoil 
are stronger as a result of reducing the numerical dissipation. However, the 
mean convergence rate with the raultigrid method deteriorates substantially. 
It is .876 with the original model and .960 with the modified (Eq. (16)) 
model. If the modification of Eq. (17) is applied, the calculated pressure 
distribution (Fig. 4a) is very close to that given in Fig. 3. Furthermore, 
the mean convergence rate of the computation is improved significantly (a 
value of .890 for 100 cycles). It should be emphasized that the function 
employed in Eq. (17) for scaling the artificial dissipation is by no means 
optimum. The lift and drag coefficients for these cases and those predicted 
with the high density mesh calculations of Ref. 7 are given in Table I. 


Table I 

Lift and drag coefficients for NACA 0012 airfoil, M = .8, a = 1.25° 

00 


Case 

Cl 

% 

Basic dissipation model 

.3330 

.0220 

Modified (Eq. 16)) dissipation 
model 

.3667 

.0235 

Modified (Eq, 17)) dissipation 
model 

.3567 

.0234 

Ref, 7 - 561 X 65 C-type mesh 

.3618 

.0236 
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Subsonlc Laminar Flow 

Numerical solutions were obtained for laminar flow past an NACA 0012 air- 
foil, The Mach number was 0.5, the Reynolds number was 5000, and the angle of 
attack was zero degrees. A C-type mesh with 256 x 64 cells (129 points on 
the airfoil) was employed in the calculations. The normal mesh spacing at the 
surface was about 6 x 10 ^ chords, and the trailing edge streamwise spacing 
was 5 X 10 chords. In Figs. 5a and 5b, the pressure and skin-friction 
distributions computed using the original dissipation model are shown. The 
absence of any pressure recovery at the trailing edge indicates the presence 
of strong viscous effects. Moreover, as denoted in Fig. 5b, the flow 
separates at the .811 chord location. There is a sudden change in the skin 
friction at the trailing edge. At least in part, this is a consequence of the 
artificial dissipation model. The convergence history for this case is pre- 
sented in Fig. 5c. In 300 multigrid cycles with the finest grid (requiring 
less than 3 minutes on the CRAY II computer) , the mean rate of convergence is 
.923. 

The surface skin-friction distribution calculated using the modified (Eq. 
(16)) model is displayed in Fig. 6a. Now, there is a significantly smaller 
decrease in the skin friction at the airfoil trailing edge. In the case of 
the modified (Eq. (17)) model, some additional scaling (with a simple second 
degree polynomial) in the streamwise direction in the immediate vicinity of 
the trailing edge was required to obtain essentially the same skin-friction 
solution and a good convergence rate. Figure 6b shows the convergence 
histories for these cases. The mean rates of convergence using the modified 
models based on Eqs. (16) and (17) are .947 and .932, respectively. 
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For all these laminar flow results the second difference dissipation 
terms are set to zero. In the shear layers, the normal physical viscous terms 
generally dominate (are an order of magnitude or more larger than) the 
numerical dissipation terms. However, even with the modified (Eq. (16)) 
model, the streamwise artificial dissipation terms are not dominated by the 
normal physical ones for a few cells surrounding a trailing edge cell. For 
this laminar case, the streamwise diffusion terms, which were neglected by the 
thin-layer approximation, may be of sufficient importance to allow domination 
of the total physical viscous effects over the artificial ones at the trailing 
edge. 

The streamlines of the recirculation zone for the modified (Eq. (16)) 
dissipation model solution are presented In Fig. 7. The longitudinal and 
lateral extents of this thin bubble are very close to those predicted with the 
other models. Figure 8 shows a velocity vector plot for this laminar flow 
problem. 

Transonic Turbulent Flow 

Nonphysical solutions for turbulent flow in the trailing edge region of 
an airfoil have been observed by many investigators. The basic factors that 
determine the accuracy of the trailing edge solution are as follows; 


1) Mesh (resolution and orientation), 

2) artificial dissipation, 

3) turbulence modelling. 
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In Ref. 33, Haase Indicates that the principal reason for inaccurate results 
is the nonalignment of the trailing-edge mesh line and streamline. Based on 
the present work and Ref. 2, this is not considered to be the main cause of 
inaccuracy. That is, qualitatively correct physical behavior can be obtained, 
even if the trailing-edge mesh line bisects the trailing-edge angle, as long 
as the artificial dissipation is sufficiently small. 

The Impact of the artificial dissipation terms is revealed in results for 

OO 

transonic flow over an RAE 2822 airfoil. In this standard test case, the 
free-stream Mach number is 0.73, the Reynolds number is 6.5 x 10^, and the 
angle of attack corrected for wind-tunnel wall effects is 2.79 degrees. The 
first set of results was computed with the basic artificial dissipation model 
and a C-type mesh having 264 x 100 cells. A view of the mesh and a blowup 
of the trailing-edge region is shown in Fig. 9. The mesh spacing in the 
normal direction at the surface is such that the first mesh point is inside 
the laminar sublayer. The spacing in the x direction at the trailing edge 
(Ax^ ) is 0.0147 chords. Figures 10a and 10b compare the pressure and 
upper surface skin-friction distributions with experimental data. Even though 
there is an adverse pressure gradient on the upper surface near the trailing 
edge, the skin friction there exhibits a substantial rise, which is not 
physically correct. 

The next set of results was obtained with the basic dissipation model and 
a finer trailing-edge mesh. The mesh for this case is presented in Fig. 11. 
The spacing Ax is 0.005 chords. This represents a reduction of almost 

a factor of three. At the shock wave, the spacing is more than twice that for 
the previous results. In Figs. 12a and 12b, the pressure and skin-friction 
variations are displayed. There is still a strong skin friction rise at the 
airfoil trailing edge. 
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The final group of results was calculated using the modified (Eq. (16)) 
dissipation model and the mesh in Fig, 11, They are shown in Figs, 13a - 
13d. The predicted pressures are in good agreement with experimental data, 
even with the coarse grid spacing in the vicinity of the shock. As indicated 
in Fig. 13b, there are two separated flow regions on the airfoil. A very 
small shock induced separation bubble occurs at about 56% chord. The trailing 
edge separation on the upper surface of the airfoil occurs approximately at 
95% chord. The behavior of the flow in the vicinity of the trailing edge is 
clearly visible in Fig. 13c. Figure 13d presents the residual and lift 
histories for the calculation. A multigrid procedure was not employed. 
Finally, in Table II the predicted lift, drag, and pitching moment 
coefficients are compared with those of experiment and Ref. 3. 


TABLE II 

Lift, drag, and pitching moment coefficients for RAE 2822 airfoil, 
M = .73, Re = 6.5 x 10^, o = 2.79° 

CO ’ 00 ' 



Cl 

S 

P 


Cd 

Cm 

Experiment (Ref. 33) 

.803 

- 

- 

.0168 

-.099 

Present (256 x 64 
C-type mesh) 

.829 

.0124 

.0051 

.0175 

-.093 

Pulliam (Ref, 3, 248 x 51 
0-type mesh) 

.824 

.0128 

.0050 

.0178 

-.092 


Cjj - pressure drag coefficient 
P 

C_ - friction drag coefficient 
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CONCLUDING REHABKS 

Improved accuracy of numerical flow solutions has been achieved by 
modifying a standard artificial dissipation model for central differencing 
schemes. With the eigenvalue scaling suggested by upwind differencing, the 
artificial dissipation in the strearawise flow direction has been reduced. 
This has resulted in a better representation of inviscid transonic flows on a 
given mesh. In addition, physically correct viscous solutions for the trail- 
ing edge of an airfoil flow have been obtained. However, the modified eigen- 
value scaling of the dissipation has resulted in slower convergence rates for 
a multigrid method driven by a multistage time stepping scheme. Improvements 
in accuracy and multigrid convergence rates have been shown possible by 
modifying the scaling with a function that depends on mesh cell aspect ratio. 
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Figure 2 
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(a) Pressure distribution. 


. Calculation of inviscid flow over NACA 0012 airfoil using basic 
artificial dissipation model (M^ = 0.8, a = 1.25 ). 
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(b) Convergence history. 
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(a) Pressure distribution. 


Figure 3. Calculation of inviscid flow over NACA 0012 airfoil using modified 
(Eq. (16)) artificial dissipation model (M^ = 0.8, a = 1.25). 
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(b) Convergence history. 
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(a) Pressure distribution. 


Figure 4. Calculation of inviscid flow over NACA 0012 airfoil using modified 
(Eq. (17)) artificial dissipation model (M^ = 0.8, a = 1.25 ). 
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(b) Convergence history. 


.5 


Figure 5 


(a) Pressure distribution. 


• Calculation of laminar flow over NACA 0012 
artificial dissipation model (M = 0.8, Re^ 


airfoil using 
■ 5000, a = 0 ] 
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lierations 


(c) Convergence history. 
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(a) Skin-friction distribution. 


Figure 6. Calculations of laminar flow over NACA 0012 airfoil using modified 
artificial dissipation models (M = 0.5, Re = 5000, a = 0^)* 
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Iterations 


(b) Convergence histories., 

A - Modified (Eq. (16)) model 
B - Modified (Eq. (17)) model 
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Figure 9a 


Partial view of mesh (264 x 100 cells) for RAE 2822 airfoil 



Figure 9b. Blowup of mesh in trailing edge region. 
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Figure 




(a) Pressure distribution. 


10. Calculation of turbulent flow over RAE 2822 airfoil with basic 

artificial dissipation model (M = 0.73, Re = 6.5 x 10^, 
a = 2.79°). « ‘ CO 


lER-STOKES 



(b) Upper surface skin-friction 
distribution. 
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Figure 11a. Partial view of mesh (256 x 64 cells) for RAE 2822 airfoil; mesh 
refinement at trailing edge. 



Figure 11b. Blowup of mesh in trailing edge region. 






(a) Pressure distribution. 


Figure 12. Calculation of turbulent flow over RAE 2822 airfoil with basic 
artificial dissipation model and mesh refinement at trailing 
edge (M = 0.73, Re = 6.5 x 10®, a = 2.79°). 
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□ EXPERIMENT, COOK 
— NAVIER-STOKES 



(b) Upper surface skin-friction 
distribution. 




(a) Pressure distribution. 


Figure 13. Calculation of turbulent flow over RAE 2822 airfoil with modified 
(Eq. 16)) artificial dissipation model ^d mesh refinement at 
trailing edge (M = 0.73, Re = 6.5 x 10®, o = 2.79°). 
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